NASA/TM-2005-213452 



Mathematical Metaphors: Problem Reformulation and 
Analysis Strategies 

David E. Thompson 
Ames Research Center 


National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California, 94035-1000 


April 2005 


The NASA STI Program Office ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA Scientific and Technical Information (STI) 
Program Office plays a key part in helping NASA 
maintain this important role. 

The NASA STI Program Office is operated by 
Langley Research Center, the Lead Center for 
NASA’s scientific and technical information. The 
NASA STI Program Office provides access to the 
NASA STI Database, the largest collection of 
aeronautical and space science STI in the world. 
The Program Office is also NASA’s institutional 
mechanism for disseminating the results of its 
research and development activities. These results 
are published by NASA in the NASA STI Report 
Series, which includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or theoreti- 
cal analysis. Includes compilations of significant 
scientific and technical data and information 
deemed to be of continuing reference value. 
NASA’s counterpart of peer-reviewed formal 
professional papers but has less stringent 
limitations on manuscript length and extent 

of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific and 
technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical confer- 
ences, symposia, seminars, or other meetings 
sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, technical, 
or historical information from NASA programs, 
projects, and missions, often concerned with 
subjects having substantial public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’ s mission. 

Specialized services that complement the STI 
Program Office’s diverse offerings include creating 
custom thesauri, building customized databases, 
organizing and publishing research results . . . even 
providing videos. 

For more information about the NASA STI 
Program Office, see the following: 

• Access the NASA STI Program Home Page at 
http://www. sti. nasa.gov 

• E-mail your question via the Internet to 
help @ sti.nasa.gov 

• Fax your question to the NASA Access Help 
Desk at (301) 621-0134 

• Telephone the NASA Access Help Desk at 
(301) 621-0390 

• Write to: 

NASA Access Help Desk 
NASA Center for AeroSpace Information 
7121 Standard Drive 
Hanover, MD 21076-1320 


NASA/TM-2005-213452 



Mathematical Metaphors: Problem Reformulation and 
Analysis Strategies 


David E. Thompson 
Ames Research Center 


National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California, 94035-1000 


April 2005 


Available from: 


NASA Center for AeroSpace Information 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301)621-0390 


National Technical Information Service 
5285 Port Royal Road 
Springfield, VA 22161 
(703) 487-4650 


Mathematical Metaphors: 

Problem Reformulations and Analysis Strategies 

David E. Thompson 

NASA Ames Research Center 
Mail Stop 269- 1 
Moffett Field, CA 94035-1000 
dethompson@mail.arc.nasa.gov 


Abstract 

This paper addresses the critical need for the development of intelligent or assisting software tools for the scientist 
who is working in the initial problem formulation and mathematical model representation stage of research. In 
particular, examples of that representation in fluid dynamics and instability theory are discussed. The creation of 
a mathematical model that is ready for application of certain solution strategies requires extensive symbolic 
manipulation of the original mathematical model. These manipulations can be as simple as term reordering or as 
complicated as discovery of various symmetry groups embodied in the equations, whereby Backlund-type 
transformations create new determining equations and integrability conditions or create differential Grobner bases 
that are then solved in place of the original nonlinear PDEs. Several examples are presented of the kinds of 
problem formulations and transforms that can be frequently encountered in model representation for fluids 
problems. The capability of intelligently automating these types of transforms, available prior to actual 
mathematical solution, is advocated. Physical meaning and assumption-understanding can then be propagated 
through the mathematical transformations, allowing for explicit strategy development. 


Introduction 

Scientific computing requires access to the sophisticated numerical manipulation and computation systems. But 
scientists also require analytic tools that support various types of problem transformations and strategic planning for 
analysis. Symbolic reasoning is excellent for management of the flow of computational components and for 
developing alternative computational strategies as might occur in a planning system; but such tools are not yet 
available for computational physics driving active models and boundary conditions. Completely apart from 
mathematical solution or simulation, both of these features of problem manipulation and symbolic reasoning should 
be brought to bear at the initial problem formulation stage. One must allow altering of mathematical 
representations, by inverting various dependency relations, by carrying out sundry symbolic substitutions and 
reformulations through differential operations, by reordering partial differential equations according to their symmetry 
bases, and by maintaining semantics throughout these transformations — only then can actual solution begin. 

There are several stages of activity involved in computational science. The first of these begins with the scientist 
finding and expressing the simplest mathematical model that is believed to capture all the detail necessary and 
sufficient to appropriately simulate or model the physical process being addressed in the problem statement. In fact, 
at a basic level, this first stage is simply defining and refining the problem statement itself. Once a model has been 
created, then solution strategies are developed which are constrained by the limits of the mathematical description, by 
the assumptions made to insure computational tractability, and by the computational resource requirements as 
available on the particular processor architectures or distributed computing environments. Various (possibly 
alternative) strategies are then selected and implemented through a suitable computational language that supports the 
algebraic, differential and difference calculus used, and that supports the symbolic manipulation of the computational 
objects themselves in order to affect linkages between various canonical computational modules. Of course, this is 
followed by debugging, by problem or theory refinement, by general exploration or sensitivity analysis of the model 
itself, and eventually by publication of some results. Of these levels of computational activity, the most basic level 
is probably the least supported by intelligent assisting software tools. Yet this most basic level of (1) finding and 
expressing the appropriate mathematical model representation, (2) analyzing or ensuring its capture of essential 
physics, and (3) transforming it to a structure amenable for applying known solution strategies, is frequently a 
highly time-consuming enterprise, fraught with error and redefinition. It could be radically shortened and verified 
through use of reliable equation-transformation operation tools. 

Most scientists carry out such initial activity off-line, and bring a more or less well-defined or reasonably- 
expressed problem statement to their workstations. However, it would be most useful if an interactive, sophisticated 


symbolic and numerical manipulation and planning system could be available to scientists at this early stage of 
research. 

This paper considers problems in fluid dynamics and instability theory. It should be noted, however, that, it does 
not refer to problems in applied computational fluid dynamics. In that domain, simulations of complex flow fields 
are created around fixed aerodynamic structures in order to observe non-laminar flow fields and to define how to 
minimize these vortex sheets by optimally redesigning the structure. These problems are structurally well-defined 
from the outset, and many years of experience have gone into creating advanced simulation and exploration tools for 
these researchers. Rather, here the focus is on the scientific understanding of the conditions under which, for 
example, various forms of turbulent transport and energy dissipation in naturally occurring fluids such as 
atmospheres or planetary interiors, might grow or transform to finite amplitude secondary flow structures; or how 
the conditions in these fluids respond to perturbations in their flow fields so that the fluid evolves into new flow 
regimes. The concentration is on the fluid physics rather than on a physical structure; the concentration is on 
maintaining understanding and insight as mathematical assumptions and transformations propagate during analysis. 

The initial representation of this type of fluid dynamics problem may be a fairly broad model description at the 
level of the relevant equations of motion, conservation laws, constitutive equations, and boundary conditions that 
embed these equations into a particular problem statement. But an actual "working" problem statement arises only 
from successive refinements, coordinate transformations, transformations in and out of various phase space 
representations, substitutions of truncated power series expansions or of new parameters, and physically or 
mathematically justified constraints or assumptions on this initial problem statement. In a sense, the general 
problem is reduced to a particular instance upon which an entire intellectual fabric may hinge. For example, later in 
this paper [and in the Appendix], a problem is outlined on isostatic recovery flow of the Earth’s mantle following 
deglaciation of the Pleistocene ice sheets under weakly non-Newtonian rheology. One starts with equations for 
conservation of momentum, of mass, and with an expression for nonlinear viscosity. After much transformation, 
substitution, and differential operation on these equations, one realizes that the essential question to be addressed is to 
clarify the actual role of coupling of harmonic disturbances in nonlinear stress-dependent rheology and induced 
wavelength harmonics. Boundary conditions are then created for that problem, and that is the ultimate problem that 
is analyzed numerically. 

The realization is an act of insight that cannot be anticipated at the outset of the problem definition. It arises 
because of the insight gained from understanding the constraints against solution methodologies embodied in the 
original problem statement, the assumptions, and selected problem reformulation strategies. It is a kind of insight 
not eagerly allocated to "assisting" systems by scientists. In particular, it is a different kind of insight than that 
gained from detailed numerical simulation and experimentation on the "whole" model, a process frequently allocated 
to high speed, vector processor architectures. It is also a different form of insight from that which arises under 
analysis of phase portraits of dynamical systems. Phase portraits are not currently tailored to expose the detailed 
physical processes responsible for transitions between bifurcated regions of the phase space. Rather, they identify 
variations in dynamical systems behavior as gleaned from identification of initial conditions of trajectories in the 
phase space, or from analysis of the evolution of eigenvalues in the complex plane as particular bifurcations 
parameters are varied with respect to characteristic flow parameters. However, the symmetries that arise from such 
analysis are the same that occur in the original differential equations, and so there should be a way of mapping 
between these kinds of representations to achieve insight on the physical process. 

Language and Metaphor 

So, why is this paper entitled "Mathematical Metaphors"? Consider a language. It has a grammar by which we 
understand the deeper linguistic structure. At the surface, the words are used to define and individuate concepts, and 
the structure allows these concepts to be interrelated. Once we have learned the grammar, and at some level the 
linguistic structure, we study the literature. The literature is a transformation of ideas through metaphor; 
understanding in the language requires metaphor because this is how the unfamiliar becomes familiar, and how new 
connections and relations are established. We use the language to reformulate the language. 

Mathematics is also a language with a surface grammar and a deep grammatical structure. We spend many years 
learning the surface grammar, which is basically manipulation of symbols for identities — everything through 
differential and integral calculus. The surface grammar is used to manipulate symbols that define and identify 
substitution possibilities wherein an equation says "the result of carrying out the calculation on the left side of the 
equals sign is the same as the result of carrying out the one on the right side." Once we get to ordinary and partial 
differential equations, to higher order and to nonlinear equations, or to variational calculus, we are at a new level of 
use in the language. There are no longer any fixed rules about what counts as a solution: a solution is any function 
that can be found to satisfy the equation, given the particular problem conditions. The bulk of scientific 
computation is this calculation by algebraic substitution of values for variables, and the substitution of calculation 
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process for equivalent calculation process. But in science theory formation or in problem statement formulation, 
prior to these substitutions of calculations for calculations, is the metaphorical transformation of equation and model 
representations, operating at the level of the structure of the language, and with its attendant insight. Our theory and 
insight is enhanced by the metaphors we create through this analytic transformation of the problem statement. This 
act is not part of a 'solution methodology' but rather is transformation of problem representation while maintaining 
physical understanding. Just as a theory is a metaphor between model and reality, so is a problem refinement a 
metaphor between the original model and a solution. This level of theory formation is the level of the expressibility 
and subtlety of the language. It is this metaphorical level of expressibility in assisting tools that is advocated. 

Discussion and Examples 

In this section, several examples are presented of the types of equation transformations that occur regularly in 
nonlinear fluids problems and in dynamical systems generally, thereby hopefully lending clarification to the ideas of 
the previous section. It is difficult to rank these operations according to which are the most useful or the most 
necessary for near-term intelligent automation. Rather, it is hoped that some of these techniques will seem within 
the grasp of some of the Intelligent Systems Research Community, or that these examples will stimulate a research 
domain for intelligent mathematical and scientific tool development. 

I. Burgers’ Equation Transformations 

First, consider Burgers' equation [Zwillinger, 1998, p.387-390], and look at two possible transformations of it that 
can make it solvable using alternative methods. The appropriate form can well depend on the physical setting in 
which the equation arises in a given problem, so one form may be more "intuitive" than another. Burgers' equation 
is generally thought of as one form of mathematical model for turbulent flow, and it is also used in the approximate 
theory for weak stationary shock waves in real fluids. It has the general form of a quasi-linear PDE, L[u] = 0, and 
with X a parameter, it has the form: 


u t + UU X - Xu xx = 0 (1) 

We can now construct what is called the Hopf transformation, in two steps (Backlund transformations generally). 
Introduce a new variable, v, and set u = v x . Substitution yields the equation: 

v x t + v x v xx - Xv xxx = 0 (2) 

and one integration by x yields: 

v t + (l/2)v x 2 - Xv xx = 0 (3) 

This equation looks worse than before, but if we then introduce a second transform of variables between v and a new 
variable w such that v = -2 X log(w) , equation (3) reduces to a simple one-dimension diffusion equation 

wt - Xw xx = 0 (4) 

which is solvable by known methods. Of course, a solution for w yields a solution for v and hence u. 

The combined transformation: 


u = -2X [log(w)] x (5) 

is called the Hopf transformation. It is presented here in two stages, because if one stops at equation (3), there is an 
alternative method of solution. Furthermore, perhaps the physical problem under consideration yields this form of 
the transformed Burgers' equation, and the diffusion solution is not what is wanted. Equation (3) can also be solved 
by so-called "separation" methods. Here, a typical solution is sought in the form: 

v = f[t + g(x)] = f(w) (6) 


for unknown f and g . Here, w = t + g(x). Upon substitution of (6) into (3), a relation arises for f and g: 



f'(w)[l - W(x)] = [g'(x)] 2 [X.f"(w) - (l/2)(f'(w)) 2 ] 


(7) 


and this can be separated into: 


[1 - Xg"(x)] / (g'(x)) 2 = [Xf"(w) - (l/2)( f’(w)) 2 ] / f'(w) 


= c [a constant] 


( 8 ) 


Hence, f and g must satisfy the system of equations: 

/.f" - (l/2)(f ) 2 - cf' = 0 
Xg" + c(g’) 2 = 1 


(9) 


These equations are integrable by reduction of order. The successes from using transformations of the form (6) for 
Burgers' equation above has led to the generalization of such transforms for use in solutions of the Navier-Stokes 
equation [see II.], whereby the stream function is assumed to have the form 

if = f[t + g(x) + h(y)] (10) 


and this then opens up another area of analysis for instability and turbulence. 


This first example with Burgers’ equation is fairly straight-forward, and is well known in the literature. Here, it 
introduces the flavor of transformation of equations tailored to specific needs. The point is that if access to suites of 
transformations were available in a computational system, then various experiments could be run by operating on 
the defining equations and discovering the symmetries that emerged. This would of course require not only the 
possibility of functional substitution, but also the availability of differential operation on equations as a whole so as 
to separate or project forms on subspaces or create, say, poloidal and toroidal component equations, and so on. 


II. Navier-Stokes Equation Transformations 

As a second example, consider some simple variable transformations on the Navier-Stokes equation. This 
equation represents conservation of momentum in fluids [Batchelor, 1967, p.147; Schlichting and Gersten, 2003, 
p.68 & lOlff.], and for incompressible fluids, it takes the vector-indices form: 


3uj/3t + uj dui/dxj = - 1/p dp/dx + v [d 2 ui/dxj3xj] (11) 

where indices i,j = 1,2,3, and uj is velocity, p is density, p is pressure, and v is constant kinematic viscosity. 
Because the fluid is incompressible, the conservation of mass is expressed by a simplified form of the continuity 
equation, namely, dufc/dxk = 0. That means that, for the two-dimensional problem (u,v; x,y), we can define a 
stream function that satisfies continuity, whereby now: 


u = r|)y and v = - r[> x ; with xpy X - xp X y = 0 (12) 

Thus, one can see that the Navier-Stokes equations could be written in terms of the stream function x[> and the 
pressure gradients by substituting the above definitions, and carrying out the necessary differentiations. Normally, if 
one is proceeding along this path, one also then cross-differentiates the equations and subtracts to eliminate pressure, 
and the resulting equation is a fourth order PDE in the stream function. The partial derivatives have essentially 
projected the original equations onto a one family set of curves in "stream-function" space. But, instead of pursuing 
that path here (and because the final example carries it out in more detail), consider the transformation that arises 
from assuming a solution of the form: 


ui = -2v/<|> [d<|>/dxi] (13) 

If this transformation is substituted into the original Navier-Stokes equation for uj , the resulting equation is a 
linear diffusion equation in (|) : 
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(14) 


3<|)/dt = v[32(|)/axi3xj] + (p(t,x 1 ,X2,X3)/2pcl)} cj) 

This equation can now be thought of as a model for viscous flow in a pure initial-value problem where pressure is 
prescribed. In this case, it turns out that the velocity field must satisfy a source distribution of 

S(t, xi, X2, X 3 ) = - 2v[3 In ^/dxjdxj] (15) 

in order to guarantee that conservation of mass is still satisfied. If, however, S = 0, then the linear diffusion 
equation for (|), (14) above, along with mass conservation, transforms into a Bernoulli equation, which classically 
relates flow to pressure and density. On the other hand, the nonlinear form of Bernoulli equation 

30/3t + 1/2 [30/3x1 30/3xj] + p/p = 0 (16) 

goes back to a linear diffusion equation by using mass conservation and the transform 0 = In <f> . 

Each of these manipulations are in the general class of reducible equations and are found in nonlinear PDE texts; 
they may provide guidance to scientists who are working on finding a reasonable representation for their particular 
problem. Such transforms and others of a similar nature would be valuable if they were available in the "working" 
knowledge of a computationally assisting software system. 

So much for these kinds of examples: where's all the insight that is supposed to arise from these transformations? 
Where's the metaphor ? To see that, one needs to consider an application of these types of transformations to real 
problems. What is outlined now in this final example is a geophysical fluid dynamics problem. Section III 
basically outlines the structure and evolution of the analysis; the actual detailed analysis, evolution, and 
transformation of this particular geophysics problem is presented in the Appendix. The problem is not solved; 
rather, the problem is transformed to a solvable form, and the appropriate boundary conditions are integrated into the 
analysis. That final analysis is then ready for programming to a numerical solution or simulation, and the exact 
physical assumptions inherent are carried forward into that simulation. 

III. Example: Isostatic Rebound of the Earth’s Mantle under Non-Newtonian Rheology 

Suppose the general problem is one of understanding the nonlinear viscosity distribution in the Earth's interior. 
There are several sources of data that can constrain such analysis: earthquake-seismic, laboratory data on phase 

transitions in mantle-like rocks, gravity anomalies over large regions of the Earth, historical records on sea-level 
change, earthquake free-oscillation data, and changes in length of day due to tidal de-spinning of the Earth-Moon 
system, to name a few. It may turn out that most of this data can be supported by more simple layered linear 
viscosity; but for now suppose one has reason to hypothesize a nonlinear viscosity distribution. Then here is an 
example strategy with transformations. 

The problem to be considered is the flow of the mantle of the Earth in response to the unloading (melting and 
removal) of the Pleistocene ice sheets about 7K to 10K years ago. We want to determine the viscosity profile in the 
mantle p(x,y) based on this flow v(x,y) of the mantle in response to the unloading stresses. Bold terms are vector 
quantities. The basic equation of motion for rebound flow in this medium is the variable viscosity form of the 
Navier-Stokes equation, representing conservation of momentum. In vector form this equation is: 

p Dv/Dt = pg - Vp + [V- (2pV)]v + V x (pV x v) (17) 

Conservation of mass is simply div v = 0, as before. The inertial terms in the convective derivative Dv/Dt can be 
ignored in this problem because they are small compared to the viscous terms and restoring stresses. It is assumed 
that the flow is confined to the (x,y) plane, with y positive vertically upward from the non-deformed surface. The 
next usual operation is to introduce a term for the vorticity as curl v into the equations of motion, which then 
become: 


V(p + 4>) = - pV x | + 2(Vp • V)v +(Vp) x | 
with % = V x v 


( 18 ) 



Here the gravity force is written as the gradient of a potential <t> . This equation can now be operated on by taking 
the curl of the whole equation (ie, V x : ) which eliminates the potential terms and projects the vector equation onto 
the (x,y) plane. This operation also serves to eliminate several of the terms in the expanded form of the differential 
vector operator relations. The vorticity vector itself points in the transverse plane, and its components lie entirely in 
(x,y). After carrying out all the operations and imposing mass conservation, the two-dimensional form of the single 
equation of motion can be written explicitly in terms of u, v, and u. with partial differentiation in x, y: 

l'( v xxx ■ u xxy + v xyy ■ u yyy) 

+ 2[p. x (v X x - u xy) + ByOxy - Uyy)] 

+ 2p. X y(Vy - u x ) + (q xx - q yy )(v x + u y ) = 0 (19) 

The horizontal and vertical velocity components can again be represented in terms of a Stokes stream function 
ip(x,y) 


u = ip y and v = - xp x (20) 

which becomes the single dependent variable describing the flow. Substitution yields: 

iK^xxxx + 2x|J X xyy + H'yyyy) 

+ 2[[X x (xp xxx + rp xyy ) + U y (l[) XXy + tpyyy)] 

+ 4pi xy ap xy + (q xx - p,yy)(lp X x - rjjyy) = 0 (21) 

Solution of this fourth-order, viscous-dominant equation of motion is sought in terms of r|j(x,y). This equation is 
equivalent to Brennen's equation for linear viscosity [Brennen, 1974, p. 3995, his eqn. (9)], but (21) includes both 
horizontal and vertical variations in the viscosity. Integrations of would then give u and v, and if an explicit 
relation is made for the viscosity in terms of the flow strain-rates, then this viscosity profile can be found as well. 
Conversely, a viscosity form could be assumed in order to start the analysis, and iterations would then converge to a 
self-consistent system. 

The difficulty is that the viscosity u.(x,y) is now a complicated function of the stress field or strain-rate field in the 
medium; hence it is also a function of r|>(x,y). The equation of motion is thus highly nonlinear in rp(x,y), and exact 
admissible solutions are not known. One cannot assume a simple modal solution for ip , such as a superposition of 
Fourier components, because in the nonlinear problem, as the deformed surface rebounds, a new stress field is 
induced by which given disturbance harmonics will couple and induce other harmonics. In addition, the change in 
strain rates during recovery will alter the viscosity profile sensed by the rebounding nonlinear fluid. That is, u must 
be considered as the functional u[i|)(x,y)]. Hence, the form of a general solution to this equation of motion must be 
derived which is appropriate for describing the rebounding flow in such a non-Newtonian medium. Transformations 
(math metaphors) become necessary .because one hopes to preserve physical understanding at the same time. Once a 
general solution if (x,y) is found, it can be used in conjunction with various boundary conditions to develop an 
equation which describes the time-change (relaxation) of the deformed surface. This rate can be compared with 
projections from data for verification of the model. 

In order to derive the general solution ip(x,y), it is necessary to impose physically reasonable assumptions on the 
nature of the viscosity variation, and to indicate the extent to which flow in the medium is coupled to this viscosity 
variation. Ideally, these assumptions should also yield mathematical tractability! One relation to select is a 
viscosity variation that can be described as weakly spatially coupled (i.e., the viscosity at one position is dependent 
on the viscosity of material nearby). One also thus has to assume that the stream function behaves smoothly over a 
narrow bandwidth about some prescribed wavelength, even though the overall variation across the spectrum may be 
significant. These assumptions are converted to mathematical constraints on the form of the solution ip(x,y) and on 
the form of the transformation analysis during reduction of the equations. It turns out that this approach provides 
both significant physical insight into the nonlinear rheology problem as well as mathematically tractable analysis. 

An acceptable method to advance beyond the statement of this fourth-order equation is to Fourier Transform it, 
solve the resulting equation in the (k,h) transform domain for W (k,h), then inverse transform this solution back to 
yield an analytic form for rp(x,y). This analytic form of ip(x,y) can then be differentiated appropriately, now with 
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some understanding of the role of coupling of harmonics, and an actual relaxation equation can be developed for the 
rate of change of the deformed surface under specified initial displacement or stress conditions. By using our r[i(x,y), 
and combining the expression for vanishing normal stress at the deformed free surface with a kinematic surface 
condition that ties the vertical velocity of that surface to the rebounding motion there, a partial differential equation 
is derived which describes the time evolution and rebound of the initially specified depressed land. Solution of this 
relaxation equation is the ultimate goal; it can be carried out numerically using finite-difference or multi-grid 
techniques. 

All work prior to articulating this relaxation equation can be considered symbolic manipulation and exploration in 
theory formation and problem statement reformulation. In this summary, the presentation has been simplified, in 
that the analysis in the Fourier Transform domain turns out to be rather involved. The convolutions integrals that 
are solved in (k,h) are of course the various products in (x,y), such as PxU'xyy, and solution of these integrals is 
what requires the stringent assumption on the form of nearest neighbor "awareness" of the viscosity relations. In 
several instances, truncated series expansions are needed, and these are substituted (more metaphor). At one point, an 
interesting first-order PDE arises to be solved for the transformed stream function W (k,h) : 

A(k,h) d W (k,h)/dk + B(k,h) d»F(k,h)/dh + C(k,h) *F(k,h) = 0 (22) 

This equation is then solved along characteristic paths in (k,h) space by a transformation of variables. The 
subsidiary equation 


dh/dk = B(k,h)/A(k,h) = function(k.h,k',h') (23) 

defines the characteristic paths in (k,h) along which the solution must exist, and this yields a constraint on the 
relation between the two wavenumbers k and h (representing the transform from x and y). The primed k and h are 
not derivatives, but rather here represent the width of the pulse for li in the convolution integrals in the k and h 
directions. This constraint dictates the form of the final solution, but it also yields insight into the nature of the 
coupling of the harmonics from the relaxing deformation. For a suitable transformation of variables, the first-order 
equation above becomes: 

A&0) dV(&e)/dZ + C(L0) f'(LO) = 0 (24) 

In this problem, the appropriate transformation of variables turns out to be; 

0(k,h) = hk' b/a and £(k,h) = k (25) 

where b and a are also functions of the pulse width. This transformation leads to a solution of the first-order 
equation in CC,0) space, and inverse transforming this yields a solution for W (k,h). Finally, the W (k,h) solution is 
inverse Fourier Transformed to yield the desired solution for our original stream function i|) (x,y). Plus, the 
nonlinearity of the viscosity is completely incorporated in this solution due to the convolution integrals. 

This i|’(x,y) can now be analytically differentiated, and finally an equation can be developed for the relaxation of the 
surface (see equation (A72) in the Appendix). It is this final equation that becomes the particular problem to be 
solved in this physical scenario. All the rest is preparation for that. 

Summary 

A wide variety of symbolic manipulation of equations can occur during the initial problem statement formulation 
stage in a theoretical applied math problem. Frequently, scientists bring fully formed problem statements to their 
workstation to seek numerical solutions or symbolically guided variation of parameters of a well-formed problem. 
However, in many of those cases, a tremendous amount of off-line analysis has been accomplished ahead of runtime, 
and this work can be referred to as problem statement formulation; starting from general equations of motion and 
transforming and reducing the equations to a solvable representation. This enterprise is in urgent need of intelligent 
mathematical manipulation tools by which the scientist can explore various transformation possibilities and see 
what physical constraints are intimately tied to each kind of transformation or assumption. This theory formation 
stage of theoretical work requires the capability to perform substitution of variables, substitution of truncated series 
representations for variables, changes in independent variables so as to reformulate equations in alternative spaces, 
carrying out integrations such as Fourier Transforms so that analytic results can be viewed and manipulated, inverse 
transforming back while retaining original meaning, and re-writing equations according to various symmetry 
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operations or from differential vector operations. In an ideal case, a user might be able to start with conservation 
laws for mass, momentum, and energy, and write them as balance equations in "english" or in pseudo-code, and later 
replace words with mathematical terms representing the processes. That would truly be metaphorical! In all of these 
operation, it is necessary to maintain a background knowledge whereby vector symbols retain meaning or body 
forces really can be considered as gradients of potentials, and so on. In addition to these intelligent data structures, 
planning systems could be exploited to manage the progress of the calculations, while retaining the experience of 
alternative strategies. Understanding would then propagate through the various transformations. Such systems 
would be very valuable and broadly used were they to come under development. 
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APPENDIX: The Nonlinear Mantle Viscosity Problem and its Transformations 


I. Derivation of the Equations of Motion and the Stream Function Solution 

This appendix addresses one particular strategy for formulating the problem of the gravitational relaxation and 
rebound of a nonlinear viscous Earth’s mantle after unloading. Insofar as one is trying to understand the physics of 
the rebound behavior, the following analytic work helps to isolate the effects of those processes that drive the global 
behavior. The fundamental exploratory work here is not to solve a rebound problem, but rather to use this rebound 
problem formulation as a domain within which to explore the impacts of various assumptions in the analysis. The 
particular strategy of the method presented here is as follows: start with expressions of conservations equations for 
mass and momentum. Conservation of mass for incompressible fluids allows the velocity field to be represented in 
terms of one variable, the stream function tf>(x,y). This stream function can also explicitly incorporate the variable 
viscosity p(x,y), because pi[xp(x,y)], and this can be inverted. The products of various derivatives of xp and p 
become convolution products in the Fourier Transformed momentum equations. The function W (k,h) that satisfies 
these equations is derived through a series of algebraic manipulations and integral transforms. Inverse transforming 
^(kjh), and using that derived t|t(x,y) in the boundary conditions allows derivation of the relaxation equation of the 
free-surface r)(x,y,t). This equation describes the time rate of change of the rebounding surface under the constraints 
and assumptions explicitly made in the particular strategy outlined here. Other strategies are possible with different 
outcomes, and these are discussed as options in the text. 

The problem to be considered is the flow of the mantle of the Earth in response to the unloading (melting) of the 
Pleistocene ice sheets about 7K to 10K years ago. The goal is to determine the viscosity profile in the mantle 
p(x,y) based on this flow v(x,y) of the mantle in response to the unloading stresses. Bold terms are vector 
quantities. The basic equation of motion for rebound flow in this medium is the variable viscosity form of the 
Navier-Stokes equation, representing conservation of momentum. In vector form this equation is: 

p[dv/dt + V • Vv] = p Dv/Dt = pg - Vp + [V • (2pV)]v + V X (pV X v) (Al) 

Conservation of mass requires div v = 0. The inertial terms in the convective derivative Dv/Dt can be ignored in 

this problem because they are small compared to the viscous terms and restoring stresses. It is assumed that the 

flow is confined to the (x,y) plane, with y positive vertically upward from the non-deformed surface. The basic 
equation of motion thus reduces to one of Stokes flow for variable viscosity: 

-V(p + pgy) -jtV X VxV + 2(V(1 • V)v + (Vp) x (V X V) = 0 (A2) 

Equation (A2) can be rewritten in terms of the vorticity, curl v, by introducing the vector 

£ = V X V (A3) 

into the equation of motion, which then becomes: 

v(p + O) = - pV x § + 2(Vp • v)v + (Vn) X § (A4) 

Here the gravity force is written as the gradient of a potential, <I>. This equation can now be operated on by taking 
the curl of the whole equation (ie, V x : ) which eliminates the potential terms and projects the vector equation 
onto the (x,y) plane. This operation also serves to eliminate several of the terms in the expanded form of the 
differential vector operator relations. Taking the curl of equation (A4) yields: 


(A5) 


0 = hV 2 § - V(Vn • 5) + 2(1 • V)V(1 + 2(V|X • V)i 
+ 2[(VjVmH)Vm] x v - §v 2 (i 
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wherein the repeated subscript m in the operator (VjVm[r)Vm is summed over spatial coordinates for each 
operator Vj . For the two-dimensional problem, all flow is confined to the (x,y) plane so that all transverse partials 
d/d z and the Eulerian velocity component w = ve 3 vanish. Furthermore, the vorticity as defined in (A3) becomes: 

| = ( dv/dx - du/dy ) e 3 = e 3 (A6) 

and points in the transverse plane so that its components lie entirely in (x,y). Here u and v are the x- and y- 
components of the velocity vector. With this restriction, the terms - V(Vp • §) and 2 • V ) V u inequation 

(A5) vanish as well. 

After carrying out all the operations and imposing mass conservation, the two-dimensional form of the single 
equation of motion can be written explicitly in terms of u, v, and li, with partial differentiation in x, y: 

^(Vxxx ■ Uxxy + Vxyy - Uyyy) 

+ 2[hx(v XX - Uxy) + Hy(v X y - U yy )] 

+ 2[i X y(Vy - U x ) + (hxx _ Hyy)( v x + Uy) = 0 (A7) 


The horizontal and vertical velocity components can be represented in terms of a Stokes stream function ip(x,y), by 
which 

U = t|ty and V = - XjJx (satisfying div v = 0). (A8) 

ip(x,y) then becomes the single dependent variable describing the incompressible flow. Substitution yields: 

xxxx + 2ipxxyy + ^yyyy) 

+ 2[hx(^ XXX + 4>xyy) + Hy(^xxy + ^yyy)] 

+ 4(t X y^xy + (Hxx - HyyM^xx - ^yy) = 0 (A9) 


This equation is equivalent to Brennen's linear viscous equation [Brennen, 1974, p. 3995, his eqn. (9)], but (A9) 
includes both horizontal and vertical variations in the viscosity. Solution of this fourth-order, viscous-dominant 
equation of motion is sought in terms of ip(x,y). Integrations of ip would then give the velocity components u and 
v, and if an explicit relation is made for the viscosity in terms of the flow strain-rates, then this viscosity profile can 
be found as well. Conversely, a viscosity form can be assumed in order to start the analysis, and iterations can 
converge to a self-consistent system. 

The problem is that the viscosity q(x,y) is now a complicated function of the stress field or strain-rate field in the 
medium; hence it is also a function of ip(x,y). The equation of motion (A9) is thus highly nonlinear in ip(x,y), and 
exact admissible solutions are not known. One cannot assume a simple modal solution for ip , such as a 
superposition of Fourier components, because in the nonlinear problem, as the deformed surface rebounds, a new 
stress field is induced by which given disturbance harmonics will couple and induce other harmonics. In addition, the 
change in strain rates during recovery will alter the viscosity profile sensed by the rebounding nonlinear fluid. That 
is, [i must be considered as the functional u[ip(x,y)]. Hence, the form of a general solution to this equation of 
motion must be derived which is appropriate for describing the rebounding flow in such a non-Newtonian medium. 
Transformations become necessary, but we hope to preserve physical understanding at the same time. Once a general 
solution ip(x,y) is found, it can be used in conjunction with various boundary conditions to develop an equation 
which describes the time-change (relaxation) of the deformed surface. This rate can be compared with projections 
from data for verification of the model. 

An acceptable method to operate on equation (A9) is to Fourier Transform it, solve the resulting equation in the 
(k,h) transform domain for the transformed stream function l I / Yk,h), then inverse transform this solution back to 


10 



yield an analytic form for rp(x,y). This analytic form of ip(x,y) can then be differentiated appropriately, now with 
some understanding of the role of coupling of harmonics in the nonlinear problem, and an actual relaxation equation 
can be developed for the rate of change of the deformed surface under specified initial displacement or stress 
conditions. By using our rp(x,y), and combining the expression for vanishing normal stress at the deformed free 
surface with a kinematic surface condition that ties the vertical velocity of that surface to the rebounding motion 
there, a partial differential equation is derived which describes the time evolution and rebound of the initially specified 
depressed land. Solution of this equation is the ultimate goal; it can be carried out numerically using finite-difference 
or multigrid techniques. 

By this analysis method, convolution integrals arise, in the transform domain, that represent the products that 
occur in equation (A9), such as pil>xxxx> MxH’xyy; and so on. These integrals, however, contain the unknown 
variable •P'lkji) for which a solution is sought. Hence, in order to evaluate the convolution integrals, either the 
behavior of y 7 ” (k,h) must be known or hypothesized in advance, or some physical assumptions must be made. 
Appropriate assumptions are discussed below. 

Let (k,h) represent the Fourier Transform variables for (x,y), and let //(k,h) and •Fik.h) represent the transformed 
variables for u(x,y) and i|i(x,y). Here we have initially assumed that li = u[i|j(x,y)] can be represented explicity 
in terms of u(x,y) . Then the Fourier transform of equation (A9) yields the following equation: 

//( k,h) * [k 4 ^(k,h) + 2k 2 h 2x F(k,h) + h 4 ^(k,h)] 

+ 2k//(k,h) * [k3y/(k,h) + kh2y/(k,h)] 

+ 2h//(k,h) * [k 2 h^(k,h) + h3«P(k,h)] 

+ 4khy6/(k,h) * kh«P(k,h) 

+ [k 2 //(k,h) - h 2 //(k,h)] * [k 2 ^(k,h) - h 2 ^(k,h)] = 0 (A10) 


Equation (A10) defines an integral equation in which the integrands are the convolved products of appropriate 
//( k . h ) and W (k,h) terms. This equation can be simplified to a solvable equation under reasonable physical 
assumptions on the nature of the viscosity variation and on the extent to which flow in the medium is coupled to 
this viscosity variation. Ideally, these assumptions should also yield mathematical tractability ! 1 have selected a 
viscosity variation that can be described as weakly spatially coupled (the viscosity at one position is dependent on 
the viscosity of material nearby), and I have also had to assume that the stream function behaves smoothly over a 
narrow bandwidth about some prescribed wavelength even though the overall variation across the spectrum may be 
significant. These assumptions are converted to mathematical constraints on the form of the solution r|)(x,y) and on 
the form of the transformation analysis during reduction of the equations. It turns out that this approach provides 
both significant physical insight into the nonlinear rheology problem as well as mathematically tractable analysis. 

If the viscosity p were constant (Newtonian rheology), its Fourier Transform would be a delta function so as to 
recover the W (k,h)-terms appropriately throughout the convolution integrals. If, instead, the viscosity IJ is 
assumed to be only slightly different from the Newtonian constant value, but is also weakly dependent at each point 
(x,y) on the viscosity at a small finite distance away, then its transform can be represented by a delta-like function: 
instead of a pure delta function at each point (k,h), //{ k . h ) would have some finite bandwidth represented by 2 k’ and 
2h' centered at each k and h. For convenience, we define the //(k,h) to be unity (normalized viscosity) in the range 
+k' and ±h' about k and h, and to be zero outside this range. This square-pulse definition approximates the 
envisioned delta-like function and acts as a square filter which when convolved with the W (k,h)-terms, greatly 
simplifies the convolution integrals. The viscosity variation is absorbed in the evaluation of the convolution 
integrals which are now non- vanishing only in a limited range. Thus, for example, we can identify 2 px'l'xyy in the 
transform domain as: 



(All) 


k+k' h+h' 

2k ji (k,h) * kh 2 W(k,h) = 2 / f (k-%)%r\ 2 W(%,r\) d^dri 

k-k' h-h' 


By restricting the half-widths k' and h' to be small enough, the transformed stream function W (k,h) can be 
assumed to vary only slightly from uniform over 2k' and 2h', whereby it may be represented as a Taylor series 
approximation in the integrand, truncated to linear terms only: 


00 00 

W (^,ri) = 2 2 (l/n!m!){d( n )/d^( n ) d( m )/dri( m ) [^(k,h)]} (§-k) n (ri-h) m 

n=0 m=0 

= W (k,h) + (^-k)dW (k,h)/d§ + (r]-h)dW (k,h)/d r] + 0(n,m>l) (A12) 


The example convolution integral above then becomes: 

2k Ji (k,h) * kh 2 *p-(k,h) = 

k+k' h+h’ 

2 / f (k-^r| 2 { W(k,h) + (^-k) dW(k,h)/d% + (r]-h) d^(k,h)/ari}d^dri (A13) 

k-k' h-h' 


By substituting the expansion (A12) into all convolution integrals, as in the above example, equation (A10) 
becomes an integral equation to be solved for W (k,h). Note that the derivatives 


dW( k,h)/31 

and 

dW(k,h)/dr\ 

are constant and integrate to 



d W (k,h)/dk 

and 

d W (k,h)/<9h 


All other functions of | and r| integrate to functions of k and h. The integrations are straightforward but tedious. 
After much algebra, the following partial differential equation results for W (k,h): 

[(8/3)(k' 3 h')(k 3 + kh 2 )] dW(k,h)/dk + [(8/3)(k'h' 3 )(h 3 + hk 2 )] a*p-(k,h)/dh 

+ [(4k'h')(k 2 + h 2 ) 2 + (4/3)(k' 3 h' - k'h' 3 )(k 2 - h 2 )] W(k,h) = 0 (A14) 


This equation can be conveniently written as: 

A(k,h) d (k,h)/dk + B(k,h) a^(k,h)/dh + C(k,h) W(k,h) = 0 


(A15) 


where: 



A(k,h) = [(8/3)(k' 3 h')(k 3 + kh 2 )] = ak(k 2 + h 2 ) 

B(k,h) = [(8/3)(k'h' 3 )(h 3 + hk 2 )] = bh(k 2 + h 2 ) 

(A16) 

C(k,h) = [(4k'h')(k 2 + h 2 ) 2 + (4/3)(k' 3 h’ - k'h' 3 )(k 2 - h 2 )] 

= c(k 2 + h 2 ) 2 + [(a-b)/2](k 2 - h 2 ) 

where a, b, and c represent the derived constant functions of k' and h'. 

Equation (A15) is then solved along characteristic paths in (k,h) space by a transformation 
subsidiary equation 

dh/dk = B(k,h)/A(k,h) = bh/ak 

defines the characteristic paths in (k,h) along which the solution of (A15) must exist. That 
defines the manner in which the solution Y (k,h) must change in k for a given change in h. 
yields the relation: 

hk-b/a = d = constant with b/a = (h'/k') 2 (A18) 

Thus, solution of the convolution integrals, such as (A12), with the above constraints from the resulting partial 
differential equation (A15), requires a stringent assumption on the form of nearest neighbor "awareness" of the 
viscosity functions and the flow stream function. That is, relation (A18) is a constraint on the relationship between 
the two wavenumbers k and h (representing the spatial transforms of x and y), and the terms k' and h'. Recall, the 
primed k and h are not derivatives, but rather represent the width of the pulse for li in the convolution integrals in 
the k and h directions. This constraint dictates the form of the final solution, but it also yields insight into the 
nature of the coupling of the harmonics from the relaxing deformation. If k' and h' were equal, then the viscosity 
nonlinearity would be isotropic in that the width of the nearest neighbor coupling would be the same in both the 
horizontal and vertical directions. Furthermore, the b/a exponent would be unity, and the wavenumbers h and k 
representing the variation of y in x and y would be within a constant of each other. Thus, we need to somehow 
determine the constant 'd' in (A18). 

Consider the constant viscosity problem to be the Newtonian limit problem of equation (A9). In that case, with 
(X = [X + E [X (1) + 0(£ 2 ) for some small parameter £, and it <(>) = constant, equation (9) reduces to: 

(ipxxxx + 2i|t X xyy + H’yyyy) = 0 (A19) 

for the Newtonian problem. Whatever the solution is to the full problem (A9), we want to be able to recover the 
solution to (A 19) in the limit of constant viscosity. This ties the relation between h and k to that which must result 
in the Newtonian limit as well. Equation (A 19) admits a modal solution of the form: 

x|t(x,y) = A exp[i(kx + hy)] (A20) 

for some arbitrary amplitude A. Substitution of this solution into (A 19) yields a dispersion relation that must hold 
between h and k 

(k 2 + h 2 ) 2 = 0 or ±(k 2 + h 2 ) = 0 or h = ±ik (A21) 

The physically reasonable solution requires that the response flow be attenuated at depth away from the deformed 
surface as e' k y , so to achieve this we must select h = +ik. Hence the particular solution of (A19) must be of the 
form 

xp(x,y) =Ae ikx e ih y =Ae ikx e- k y (A22) 


of variables. The 


(A17) 

is, equation (A17) 
Solution of (A 17) 
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But in the Newtonian limit, h' and k' also vanish uniformly, so b/a approaches unity so that the constant d in (A18) 
must equal +i. For the non-Newtonian problem, (A18) must hold for varying h and k and for varying h' and k' 
(with h'/k' not necessarily unity), but the product hk'°' a = i must always hold along characteristic paths in (k,h) 
space. This constraint is carried forward into the form of the general solution we are deriving here for rp(x,y) as a 
solution to equation (A9). 

For a certain transformation of variables, the first-order equation (A15) with (A16), that arose from solutions of 
the convolution integrals, can be transformed to a first-order equation of the form: 


A(£,0) dW(Q,Q)/dti + C(£,0) V(&Q) = 0 (A23) 


In this problem, the appropriate change of variables needed to achieve this transformation from (A15), with 
characteristics as defined in (A 18), turns out to be: 

£(k,h) = k and 0(k,h) = hk' b/a 


(A24) 


k = ^ and h = 0k+ b / a = 0^ a 

(See, for example, discussion in Dennemeyer, 1968, p. 12-15). The coefficients A(k,h) and C(k,h) in (A15) or (A16) 
become, under this transformation: 


AO) = a'Cf'C 2 + 02^2 b/a] 

CO) = cfC 2 + 02'c2b/a] + [(a-b)/2][C 2 -0 2 ^2b/a] 


(A25) 


If equation (A23) is integrated with respect to £, while holding 0 fixed, then W fC,0) must satisfy the general form: 

«FO) = f[0] exp [ - JcO) / A(C,0) d£, ] = f[0] x(^,0) (A26) 

where f[0] is an arbitrary function. The exponential integral part '/(’CO) satisfies (A23) alone, and inverse 
transforming its variables using (A24) yields a particular solution, s(k,h), of equation (A15), whereby 

s(k,h) = xK(k,h), 0(k,h)] = X [k, 0(k,h)] (A27) 

The general solution of equation (A15) is then: 

W (k,h) = f[0(k,h)] s(k,h) (A28) 


where f[0] is arbitrary. 

Accordingly, solution of equation (A 15) is reduced to solving equation (A26) for the exponential integral part 

Xfe0) = exp [ -fc(t,d) / A(£,0) d£ ] , (A29) 

substituting the result into equation (A27), and re-transforming the variables back from (£,,0) to (k,h) using the 
relations (A24). A form of f[0] is then selected, and that function times s(k,h) yields the general solution desired. 
Continuing the strategy, when the f^k.h) solution (A28) is inverse Fourier Transformed, it yields the desired 
solution for our original stream function x)>(x,y), solving equation (A9). Plus, the nonlinearity of the viscosity is 
completely incorporated in this solution due to the form of the convolution integrals selected. This r)>(x,y) can then 
be analytically differentiated, and finally an equation can be developed for the relaxation of the deformed land surface. 
It is this final equation that becomes the particular problem to be solved in this physical scenario. 
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To solve (A26), equation (A29) must be solved. From (A25), the ratio of C/A becomes the needed integrand in 
equation (A29) as follows: 


Cfc0)/A&0) = |c/a|['C 2 + 0 2 (;2b/a] /£ 

+ [(a-b)/ 2 a]fc 2 -e 2 52 b/a ]/^2 + 0 2^2b/a] (A30) 


or 


-fc(Z,6) / A(^,0) = - [c/a] ft, [1 + 0 2 ^ ] 

- [(a-b)/2a] J{[1 - 0 2 ^] /^[l + 0 2 ^]} d £ (A31) 

where r = (2b/a) - 2 for notational simplification. Note that if b = a (i.e., h’ = ±k') as in the Newtonian problem, 
then r vanishes. Evaluation of the integrals in (A31) requires evaluation of four integrals: 

(i) - [c/a] dt; 

(ii) - [c/a] 0 2 f C r+I d'C 

(hi) - [(a-b)/2a] f 5' 1 /[I + 0 2 ? ] d^ 

(iv) - [(a-b)/2a] 0 2 f 'O'-' /[ 1 + 0 2 C r ] d'C 

Integrals (i), (ii), and (iv) are trivial and yield: 


(i) 

- [c/a] f 'C d'C = - [c/2a] 'C 2 


(ii) 

- [c/a] 0 2 f C'' +l d'C = - [c/2b] 0 2 C r+2 

(A3 3) 

(iv) 

- [(a-b)/2a] 0 2 f ?-l /[l + Q 2 ^ ] d^ = - (1/4) In [1 + Q 2 Q ] 



Integral (iii) is not quite so straightforward. In the case that h 1 ~ k 1 , b/a = (h’/k') 2 is close to but not equal to 

unity, and r = (2b/a) - 2 is very small but not zero. Then for small r, the function l) in integral (iii) can be 
expanded as 


(A3 2) 

with r = (2b/a) - 2 


= 1 + r ln^ + (rln^) 2 /2! + ... 


(A34) 



which for small r can be truncated to 

« 1 + rln£ (A3 5) 


Under this restriction of small r and of h' ~ k', integral (iii) in (A32) becomes approximately: 

- - [(a-b)/2a] f (U 1 /[(l + 0 2 ) + (0 2 r) In £] d£ = - [(a-b)/2a] JV 1 /[p + q In £] (A3 6) 


By letting U = In z so that dU = dz/z , and with p and q defined implicitly as in (A36), integral (iii) becomes: 

- [(a-b)/2a] J" [p + qU]' 1 dU = - [(a-b)/2a] (1/q) In (p + qU) 

or 

- [(a-b)/2a] f t,' 1 /[I + ] d^ = - [(a-b)/2a] (l/0 2 r) In (l + 0 2 + 0 2 r In £) (A3 7) 


The sum of equations (A33) and (A37) are the required solutions to the integral (A31). But the solution (A37) 
holds for small r only, and it is singular if r vanishes. Thus h' ~ k' but h' * Id. Hence, under the analytic 
formulation here, the "windows" of nearest-neighbor interaction for the viscosity, as explicitly represented in the 
convolution integrals, must be slightly different in the horizontal and vertical dimensions. The viscosity is not 
isotropic, and coupling interactions are different vertically from horizontally. 

Combining equations (A33) and (A37) yields: 

-Jcfe0) / A(£,0) d£, - - [c/2a] ^ 2 - [c/2b] 0 2 ^+ 2 - (1/4) In [1 + Q 2 Q ] 

- [(a-b)/2a] (l/0 2 r) In (l + 0 2 + 0 2 r In t) 

then using again » 1 + r In'C and the definition for r, the integral (A3 1) simplifies to: 

-fc(&0) / A(£,0) d^ « - [c/2a]^ 2 - [c/2b]0 2 ^ 2b/a + [(1 - 0 2 )/40 2 ] In [1 + 0 2 ^ r ] (A3 8) 


Working backwards, this solution can now be substituted into equation (A29) for /fC.O) yielding: 

[(1 - 0 2 )/40 2 ] 

x(?,0) = exp [ - [c/2a]^ 2 - [c/2b]0 2 ^ 2b/a ] (l + O 2 ^) (A39) 


Recalling the variable changes given in equations (A24), equation (A39) can be rewritten in terms of the original 
variables (k,h), and then it can be substituted into equation (All) for the particular solution s(k.h): 
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1(1 - ( hk' b/a ) 2 )/(2hk' b/a ) 2 ] 


s(k,h) = exp [ - [c/2a]k 2 - [c/2b]h 2 ] (l + (h 2 /k 2 )) 


(A40) 


However, because solution of the subsidiary equation (A 17) dictates that relation (A18) must hold, and that the 
product hk“ b ' /;i in this relation must be +i, the exponent in equation (A40) reduces to -1/2, yielding the particular 
solution to be: 


- 1/2 

s(k,h) = exp [ - [c/2a]k 2 - [c/2b]h 2 ] (l + (h 2 /k 2 )) (A41) 


from which is derived the general solution to equation (A15): 

- 1/2 

W (k,h) = f[0(k,h)] s(k,h) = f [hk- b/a ] exp [-[c/2a]k 2 -[c/2b]h 2 ](l + (h 2 /k 2 )) (A42) 


Recall that f [hk' b ^ a ] is arbitrary; one may select it as follows. Because the characteristics requires hk“ b,/a be 
constant along paths in (k,h), and because the boundary conditions in the Newtonian case require h = ik, so that the 
constant must be +i, select f to be the particular case of unity by choosing f as: 

f [hk-b/a] = [hk-b/a] 4 = ( +i )4 = i 

Then the general solution reduces to: 


- 1/2 

«P"(k,h) = exp [-[c/2a]k 2 -[c/2b]h 2 ](l + (h 2 /k 2 )) 


(A43) 


This is the solution to equation (A15), derived from equation (A10), which must now be inverse Fourier 
Transformed to yield a general solution to equation (A9) for rp(x,y). This is accomplished by integrating W from 
(A43) as follows: 


xjt(x,y) 


(1/2jc) 2 f f W (k,h) e ih y e ikx dh dk 


(A44) 


Now, instead of carrying out this integration immediately, consider some physical constraints in order to 
experiment with the behavior of this solution in terms of the coupling of response harmonics during rebound. 
Assume first that rp(x,y) could actually be well represented in the horizontal x-direction by superposition of 
independent Fourier components. In particular, consider the assumption in which a modal solution of the form e* kx 
in the horizontal direction is sufficient to describe the flow. Then, even though the stress field is changing as the 
rebound occurs, that change in the horizontal direction is slight compared to that changing in the vertical y-direction 
so that there is very little coupling of changing horizontal viscosity distribution to changing horizontal flow 
velocity u. This assumption leads to a simplification of integral (A44) in that we can now assume no integration 
over k, and that (perhaps within a factor of 1/2 ji) (A44) reduces to 

x|t(x,y) = G(k,y)e lkx (A45) 
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Then the half-inverse transform G(k,y) of W (k,h) , now for each fixed wavenumber k, is defined by: 

oo 

G(k,y) = (l/2it)/ W (k,h) e ih y dh 

-00 

or 

00 

G(k,y) = (1/2 tc) exp [-[c/2a]k 2 ] J'fl + (h 2 /k 2 )]" 1/2 exp [-[c/2b]h 2 ] e ih y dh (A46) 

-00 

To evaluate this integral, first expand the 1 + (h 2 /k 2 ) term as the series 

[1 + (h 2 /k 2 )]-!/ 2 = 1 - (1/2) (h/k) 2 + (1-3/2-4) (h/k) 4 - (1-3-5/2-4-6) (h/k) 6 + ••• 

or 

oo v- 1 

[1 + (h 2 /k 2 )]- 1/2 =1+2 (-l) v {n [(l+2m)/(2+2m)]} (h/k) 2v (A47) 

v= 1 m=0 

Let A(k) be defined as the function of k outside the integral over dh in (A46): 

A(k) = (l/2jt) exp [-[c/2a]k 2 ] (A48) 

Then substituting this definition and the series expansion (A47) into the integral (A46) yields the following integral: 

oo 

G(k,y)/A(k) = J'exp [-[c/2b]h 2 + ihy] dh + 

-00 

oo v- 1 00 

+ 2 (-l) v { II [( 1 +2m)/(2+2m)] } ( l/k) 2v f h 2v exp [-[c/2b]h 2 + ihy] dh (A49) 

v= 1 m=0 - 00 


The first integral of (A49) is evaluated using a Gaussian identity, whereby, let: 


oo 

J* exp [-[c/2b]h 2 + ihy] dh 

-00 


00 

: J* exp [-ah 2 + ihy] dh 

-00 


a = c/2b 


Then recognizing that 


oo 

J* e' h - v dh 

-00 


00 

J ' cos (hy) dh 

- 00 
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and that 


oo 

J* exp [-ah 2 ] cos (hy) dh = (n/a) 1/2 exp [-y 2 /4a] ; 1/a = 2b/c ; l/4a = b/2c 

-00 

By direct substitution, the first integral of (A49) becomes: 

oo 

f exp [-[c/2b]h 2 + ihy] dh = [2jrb/c] 1/2 exp [- (b/2c)y 2 ] (A50) 

-00 


The second integral of (A49) is evaluated as follows. Let p = c/2b and 2q = iy . Then it takes the form: 
00 00 
J h 2v exp [-[c/2b]h 2 + ihy] dh = J* h 2v exp [-ph 2 + 2qh] dh 


which admits the series solution (see, for example, Grobner and Hofreiter, 1966, vol 2, form (6), p. 65): 

v 

= (2v)! exp(q 2 /p) [jt/p] 1/2 (q/p) 2v 2 [l/(2v-2r)! r!] (p/4q 2 ) r ; p>0 

r=0 


Substituting back for p = c/2b and q = iy/2 gives the desired function of y as: 


oo 

J ' h 2v exp [-[c/2b]h 2 + ihy] dh = 

-00 


= [2jtb/c] 1/2 exp [- (b/2c)y 2 ] [2biy/2c] 2v 2 [(2v)!/(2v-2r)! r!] (-c4/2b4y 2 ) r 

r=0 


= [2jtb/c] 1/2 exp [- (b/2c)y 2 ] (i) 2v (b/c) 2v 2 [(2v)!/(2v-2r)! r!] (-l/2) r (b/c) _r y 2v_2r 

r=0 


= [2jtb/c] 1/2 exp [- (b/2c)y 2 ] (i) 2v 2 [(2v)!/(2v-2r)! r!] (-l/2) r (b/c) 2v ' r y 2v ' 2r (A51) 

r=0 


Then substituting solutions (A51) and (A50) appropriately into the integral solution (A49), the inverse transform of 
G(k,y), for fixed k, becomes: 


G(k,y) = A(k) [2jtb/c] 1/2 exp 


[- (b/2c)y 2 ] { 1 + I(k,y)} 


(A52) 
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where, noting that (-l) v (i)- v = 1, the series I(k,y) is: 


oo y- 1 V 

I(k,y) = 2{ II [(l+2m)/(2+2m)]}k' 2v 2 [(2v)!/(2v-2r)! r!] (-1/2 ) r (b/c) 2v - r y 2v - 2r (A53) 

v=l m=0 r=0 


and where A(k) is defined in equation (A48). Note that this representation for G(k,y) has not been transformed from 
k-space to x-space. Rather a modal dependence in the x-direction has been assumed, yielding equation (A45). The 
assumption that if (x,y) = G(k,y)e'^ x rather than the integral of G(k,y)e'^ x implies that a modal dependence of i|) 
in the x-direction is acceptable, and the complete Fourier inverse transform is not necessary. Physically, this means 
one assumes either that no new harmonics are induced in the x-direction; or that any given change in the vertical 
stress field does induce changes in the x-direction, but that the change is small enough that once these harmonics are 
defined, they can be superposed to yield the total stream function if(x,y) through iteration over the k wavenumbers. 

Note that because the original motivation for this work was to understand the actual role of coupling of the 
nonlinear, stress-dependent viscosity field to changes in the rebound response flow, and to understand the changes of 
the viscosity field due to these changes in flow stresses, it is rather artificial to now make simplifying assumptions 
such as (A45) that partially eliminates the physics that was argued as fundamental to this problem! However, in an 
effort to keep track of the effect of various approximations and their relation to the resulting physical solution, it is 
worthwhile imposing and relaxing these assumptions sequentially so that the role each plays in the response flow 
can be isolated. This procedure then becomes a form of “sensitivity analysis” on the structure of the family of 
solutions to the response flow problem. So, for example, another reasonable action might now be to return to the 
original inverse Fourier Transform, as outlined in equation (A44), and carry out the full transform rather than making 
the simplifying assumption (A45) of modal behavior in the horizontal direction. But we have not yet seen the effect 
of this assumption on a resulting relaxation equation. Hence, before assumption (A45) is replaced, it is best to 
continue in this mode so that there are final results to compare when the more complicated problem is solved. In 
addition, it is not yet clear what we wish to intend as the essence of nonlinear rheology in this problem. At a 
fundamental level, the problem involves recognizing that the viscosity in the planetary interior is assumed to be a 
function of the changing stress field; but that functional relation is yet to be specified. If the viscosity were an 
explicit function of the stress, then in order to develop the viscosity profiles, one would need to calculate the stress 
field for the deformed medium either by alternative methods, or by assuming that the flow field derived here could be 
written in terms of flow strain-rates, and that the viscosity is a function of these derived strain-rates. The latter is the 
preferable strategy, but even given that, there are still options. For example, we need to determine the specific 
functional relation for u(x,y) as a function of derivatives of i|>(x,y). Ideally, the form selected could be justified on 
the basis of either experimental rock-deformation data or on convective instability models or residual gravity 
anomalies. Alternatively, a proxy representation that mirrors the structure of the viscosity profile could be assumed, 
such as an exponential approximation. In this analysis one had to make assumptions on the nature of the coupling 
between the k and h wavenumbers, and on the square-wave approximation to nearest-neighbor coupling over the k' 
and h' intervals just so that an analytic solution of the convolution integrals could proceed. This constraint must be 
maintained in order to maintain applicability of the solutions already derived. Hence, the square wave approximation 
in the transform space implies that the viscosity must look like a product of sine functions in (x,y)-space, using 
wavenumbers k' and h'. But even that restriction does not yet yield a closed problem. It is still open to selection, 
specifically, how one wishes to indicate the coupling between the viscosity changes in the vertical flow and those in 
the horizontal flow. For example, should the nonlinearity in changes in wavenumber h induce additional vertical 
response modes or should it induce both vertical and horizontal modes, seen as changes in k as well. Similarly, 
should changes in k harmonics induce changes in k, changes in k and h, or induce no new modes at all, as is 
assumed in equation (A45). Because of all of these various complexities, the problem progresses sequentially and by 
explicit tracking of the effect of the assumptions. All that one is committed to is the belief that as a certain surface 
displacement begins to rebound, the stress field in the mantle material, to depths on the order of the scale of the 
horizontal scale of the displacement, will change. And, thus the viscosity field, which is a function of this changing 
stress field will also change, so that the form and rate of the rebound will change. This stress or strain-rate 
dependence of viscosity is what is meant by the nonlinear rheology in this problem. The evaluation of specifically 
how that dependence induces other response characteristics is the core investigation of this analysis. 
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Given all these caveats, consider again the equation derived from making assumption (A45), namely equation 
(A52) above. Use this solution for ip(x,y) to derive the time-change (relaxation) equation for the rebounding surface. 
That equation must then be solved in conjunction with boundary conditions, which are now derived. 


n. Derivation of the Boundary Conditions 

The boundary conditions to be imposed at the rebounding surface, as represented by the solution (A52) and (A53) 
for (A45) above, are: a dynamic free surface condition that determines the driving normal stress, and two kinematic 
or matching conditions. The kinematic conditions at the free surface require that the horizontal velocity vanish at the 
displace surface, and that the vertical velocity exactly correspond to the time rate of change of that displaced free 
surface, say y = q(x,t). Because horizontal flows are negligible compared with the vertical rebound, and because the 
initial vertical displacement r| is small compared to the horizontal scale of the rebounding feature, the kinematic 
conditions can be linearized to order r| and be applied at y = 0. Hence to order r|, the kinematic conditions become: 

u = dxp/dy = 0 aty = 0 (A54) 

v = -3i|t/dx = 3r|/dt aty = 0 (A55) 

Evaluation of the derivatives of i|) in both of these conditions, as well as in the normal stress condition presented 
below, requires evaluation of the derivatives of G(k,y) and hence of l(k,y) from equations (A52) and (A53). For 
convenience, the various derivatives that will be needed are presented together here. From (A52): 

dG(k,y)/dy = - A(k) [2jtb/c] 1/2 exp [- (b/2c)y 2 ] ^(by/c)[l + I(k,y)] - dl(k,y)/3yj“ (A56a) 

and 

d 2 G(k,y)/dy 2 = - A(k) [2jtb/c] 1/2 exp [- (b/2c)y 2 ] 

• < ^(b/c)(l - by 2 /c)[l + I(k,y)] + (2by/c) dl(k,y)/dy - d 2 I(k,y)/dy 2 (A56b) 

with A(k) defined in (A48); and from (A53): 

oo v- 1 

dl(k,y)/dy = 2{ II [(l+2m)/(2+2m)]}k' 2v 

v=l m=0 

v 

• 2 [(2v)!/(2v-2r)! r!] (-l/2) r (b/c) 2v - r (2v-2r)y 2v - 2r - 1 (A57a) 

r=0 


and 


oo v- 1 

d 2 I(k,y)/dy 2 = 2{ II [(l+2m)/(2+2m)]} > k‘ 2v 

v=l m=0 


• 2 [(2v)!/(2v-2r)! r!] (-l/2) r (b/c) 2v - r (2v-2r)(2v-2r-l)y 2 v-2r-2 (A57b) 

r=0 
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Using relations (A52), (A53) for our general solution (A45), and substituting in the required derivative (A56a), 
kinematic condition (A54) requires that: 

u|y=o = [dH>/dy]| y =o = e ikx [dG(k,y)/3y]| y=0 = 0 


or that 

- A(k) e lkx [2jrb/c] 1/2 exp [- (b/2c)y 2 ]| y=0 ^(by/c)[l + I(k,y)] - dl(k,y)/3y| ] y=0 = 0 (A58) 


But this relation is automatically satisfied because both 

(by/c)[l + I(k,y)] and <9I(k,y)/3y 

vanish at y = 0, for all values of v and r in the appropriate series summations for 3I(k,y)/3y. 

Condition (A55) requires that at y = 0: 

v = 3r|/3t = -3x|t/3x = -ike ikx G(k,y)| y= o (A59) 

Here, G(k,y) evaluated at y = 0 becomes simply a function of k, but it is not just the amplitude A(k) times 
[2jtb/c]l/ 2 as might be expected from first glance at equation (A52) and (A53). The series I(k,y) in (A53) does not 
vanish at y = 0 because before I(k,y) is evaluated at y = 0, the series must be expanded. When the expansion is 
carried out, there are terms in which y® appear, specifically whenever 2v = 2r. These terms are of course equal to 
unity and are no longer a function of y. They hence do not vanish when the series is evaluated for y = 0. Thus, 
from (A52): 


G(k,y)| y=0 = A(k) [2rcb/c] 1/2 exp [- (b/2c)y 2 ]| y=0 { 1 + I(k,y)| y=0 } 
= A(k) [2jtb/c] 1/2 { 1 + I(k,y)| y =o} 


(A60) 


and 


I(k,y)| y =o = - (l/2k 2 ) (b/c) + (9/8k 4 ) (b/c) 2 - (225/48k 6 ) (b/c) 3 + - 


co v- 1 

= swill [( 1 +2m) 2 /(2+2m)] } ( 1 /k) 2v (b/c) v (A61) 

v=l m=0 


Thus, with A(k) defined in (A48), the second kinematic condition (A55) reduces to: 


3r|/3t = -ike ikx A(k) [2itb/c] 1/2 

co v- 1 

• { 1 + Z(-l)v{lI [(l+2m) 2 /(2+2m)]}(l/k) 2v (b/c) v } (A62) 

v=l m=0 
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In addition to these two kinematic conditions, (A58) and (A62), the normal stress must vanish at the free, 
displaced surface y = r|(x,t). In this problem, the normal stress o yy at the surface is just the hydrostatic pressure 
there plus stresses due to the rebounding viscous flow. Hence, this condition says that: 

a yy = p(y) = 0 ony = rj (A63) 

An equation for p(y) must be determined from the vertical component of the original Navier-Stokes equation, (A2), 
using the stream function definition (A8) for the velocity field. This equation is: 

Py(y) = ~Pg " M-C^xxx + ^Pxyy) ~ PxC'Pyy ~ ^Pxx) ~ 2p,yl|t X y (A64) 

Integration over y of equation (A64) provides the required relation for p(y) which must then vanish on y = r). In 
order to integrate (A64), we must have an explicit relation for the viscosity p(x,y) and its derivatives. Because the 
viscosity has been defined in the transform domain as a square- wave filter of half- widths k’ and h’, centered on k and 
h, so that the convolution integrals could be evaluated, the viscosity in (x,y) space must be represented by a product 
of sine functions that would yield this square filter in (k,h) space. Hence, for consistency in this analysis, the 
viscosity must be represented as: 


H(x,y) = sin 2jtk ' x sin 2jlh 'y = sinc(oix) sinc(py) (A65) 

2nk'x 2ith'y 


where the wavenumbers a and [1 are a = 2nk' and (3 = 2jrh'. With this definition, the required derivatives of 
viscosity in equation (A64) become: 


u x (x,y) = sinc(|3y) {(cos ax)/x - (sin ax)/ax 2 } 
u y (x,y) = sinc(ax) {(cos |3y)/y - (sin |3y)/|3y 2 } 


(A66) 


The stream function terms in equation (A64) are defined in terms of G(k,y), whereby: 


M’xxx = -ik 3 e ikx G(k,y) 

Xjtxyy = ik e lkx d 2 G(k,y)/dy 2 
xjtyy = e lkx <9 2 G(k,y)/5y 2 
ijt xx = -k 2 e ikx G(k,y) 

x|t xy = ik e lkx 3G(k,y)/dy (A67) 

and G(k,y) is defined in (A52) and (A53), A(k) defined in (A48), and with the derivatives being presented in (A56a,b) 
and (A57a,b). Using all the relations (A67), (A66), and (A65), the integration of equation (A64) yields p(y) needed 
for the normal stress condition (A63). This integration can be written as follows: 

y y y 

p(y) = -pgy + e lkx A(k) [2jtb/c] 1/2 ^DlJ Fldy + D2J F2dy + D3J F3dy (A68) 

000 
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where the factors D and integrands F are: 


Dj = Di(x; k, k 1 , h') = ik 3 (sin ax)/ ax + k 2 (cos ax)/x - k 2 (sin ax)/ ax 2 

= ik 3 sinc(ax) + k 2 /x [(cos ax) - sinc(ax)] 

D 2 = D 2 (x; k, k', h') = ik(sin ax)/ ax - (cos ax)/x + (sin ax)/ ax 2 

= ik sinc(ax) - (l/x)[cos ax - sinc(ax)] 

D 3 = D 3 (x; k, k 1 ) = 2ik(sin ax)/ ax = 2ik sinc(ax) 


Fi = Fi(y;k, h') = sinc((3y) exp [- (b/2c)y 2 ] *£l+I(k,y)} 


(A69) 


F 2 = F 2 (y; k, h') 


sinc(|3y) exp [- (b/2c)y 2 ] 

•{(b/c)(l - by 2 /c)[l + I(k,y)] 

+ ( 2 by/c) 5I(k,y)/dy - 5 2 I(k,y)/dy 2 } 


F 3 = F 3 (y; k, h') = (l/y)[cos (3y - sinc((3y)] exp [- (b/2c)y 2 ] 

' { (by/c) [ 1 + I(k,y)] - dl(k,y)/dy}> 


Using these definitions in (A68), condition (A63) yields an equation for the displaced surface rp 

pgr) = e lkx A(k) [2itb/c] 1/2 ^Di(x; k, k', h ')f Fi(y; k, h')dy 

0 


+ 


y y 

D 2 (x; k, k, h') J F 2 (y ; k, h')dy + D 3 (x; k, k') J F 3 (y ; k, h')dy } | 


y=n 


(A70) 


or, by representing the integrals in (A70) by a general function T(x, y; k, k', h 1 ), evaluated at y = q: 

Pgr| = T(x, y; k, k', h’)| y=r , e ikx A(k) [2jrb/c] 1/2 


(A71) 
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HI. The Relaxation-Equation for Rebound under Prescribed Nonlinear Viscosity 

By inverting (A71) for the normal stress boundary condition and solving for pgr)/T , equation (A71) can be 
combined with (A62) by eliminating the common factor e*^ x A(k) [2jtb/c] ' ^ 2 . This then yields an equation for the 
time rate of change of the displaced surface: 

dq/dt = -ikpg [r] / T(x, y; k, k', h')| y= J ■ 

oo v- 1 

•{ 1 + 2(-l)v{n [(l+2m) 2 /(2+2m)]}(l/k) 2v (b/c) v } (A72) 

v=l m=0 


Equation (A72) is the relaxation equation which must be solved for r|(x,t), given an initial r|(x,0), at each lateral 
point x for a given wavenumber k. At the same time, k must be iterated before composition for a total 
wavenumber distribution, and the integrals contained in the function T(x, y; k, k', h') must be solved during each 
iteration on k and change of parameters k' and h'. 

(A72) is the equation that has been sought. It describes the time rate of change of the rebounding surface under 
the constraints and assumptions explicitly made in the particular strategy outlined here. Were the problem to have 
been one of linear viscosity, this same relaxation equation would have had the simpler form of: 

dr|/3t = [pg/2pk]r| (A73) 

This relaxation equation can be solved analytically to yield: 

q = 11 o CXp[t/x R ] . 


where the relaxation time Tr = 2llk/pg and X K = (r|) | (}T\ldi \ 

Ultimately, numerical, iterative solutions to equation (A72) could be carried out. These would yield relaxation 
times that would describe the rebound behavior and its variation over time due to the nonlinear coupling of the 
stresses with the viscosity. 
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